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Matrix determinants play an important role in data analysis, in particular when Gaussian pro¬ 
cesses are involved. Due to currently exploding data volumes, linear operations - matrices - acting 
on the data are often not accessible directly but are only represented indirectly in form of a com¬ 
puter routine. Such a routine implements the transformation a data vector undergoes under matrix 
multiplication. While efficient probing routines to estimate a matrix’s diagonal or trace, based 
solely on such computationally affordable matrix-vector multiplications, are well known and fre¬ 
quently used in signal inference, there is no stochastic estimate for its determinant. We introduce 
a probing method for the logarithm of a determinant of a linear operator. Our method rests upon 
a reformulation of the log-determinant by an integral representation and the transformation of the 
involved terms into stochastic expressions. This stochastic determinant determination enables large- 
size applications in Bayesian inference, in particular evidence calculations, model comparison, and 
posterior determination. 
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I. MOTIVATION 


Current and future physical observations generate huge 
data streams to be analyzed. Particle physics, biophysics, 
astronomy, and cosmology are representatives of current 
scientific fields of interest that are undergoing a revo¬ 
lution driven by increasing data volume. Typical large 
data sets in cosmology are, for instance, observations of 
the cosmic microwave background mm as well as of the 
large-scale structure [Si] as they are often wide- or all¬ 
sky observations carried out by telescopes with remark¬ 
able resolution. In order to extract information about the 
universe or physics in general, Bayesian inference meth¬ 
ods becomes more and more frequently used as their large 
computational demands become more feasible thanks to 
technology developments. The signal of interest to be 
extracted from data could be almost everything, ranging 
from just a single parameter (e.g., the level of local non- 
Gaussianity of the cosmic microwave background mm) 
to a full four-dimensional reconstruction of the structure 
growth in the universe 0191. Such ambitious Bayesian 
analyses often invoke linear transformations of the data 
or of estimated signal vectors. 

The size of the involved data and signal spaces of¬ 
ten bans the explicit representation of matrices acting 
on these spaces by their individual matrix elements. A 
prominent example appearing in many analyses is, for in¬ 
stance, the covariance matrix of a multivariate Gaussian 
distribution of a vector valued quantity, which describes 
the two-point correlation structure of the said quantity. 
Due to their large dimensions such matrices are often 
only representable by a computer routine, which imple¬ 
ments the application of the matrix to a vector with¬ 
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out storing or even calculating the individual matrix el¬ 
ements. Such routines often invoke fast Fourier transfor¬ 
mations and other efficient operations, which in combi¬ 
nation render nonsparse matrices into easily computable 
basis systems. We refer to such a matrix as an implicit 
matrix. For instance, calculating the model evidence of¬ 
ten requires calculating determinants of such matrices. 
This work provides an efficient way to numerically calcu¬ 
late determinants given only by an implicit matrix rep¬ 
resentation. 

The remainder of this work is organized as follows: In 
Sec. nil we introduce the formalism of the stochastic esti¬ 
mation of an implicit matrix and present two numerical 
examples. Section |III| provides a perspective of possi¬ 
ble applications in science. Results are summarized in 
Sec.E 


II. PROBING THE LOG-DETERMINANT OF 
AN IMPLICIT MATRIX 

A. Formalism 


Let A = (ttij) € be an implicitly defined, 

complex-valued, square matrix of order n. Implicitly 
means that the particular entries of the matrix are not 
accessible, for instance, if dealing with large data sets, 
where an explicit storage of A might exceed the memory 
of the computer. However, the action of the matrix as a 
linear operator is assumed to be known and given by a 
computer routine implementing the mapping x Ax. 

Motivated by applications in science and statistics 
(Secs. |l] and III), in particular by signal reconstruction 
techniques and model comparison in astronomy and cos¬ 
mology, where the determinant of a covariance matrix 
is required (Sec. Ill), we constrain the variety of differ¬ 


ent types of matrices by requesting that the matrix A of 
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interest is either weak diagonal dominant or Hermitian 
positive definite. The term weak diagonal dominant is 
defined by 

^ ^ ( 1 ) 

while Hermitian positive definite means 

= A and x^Ax >0 Vx G C"\{0} (2) 


(spectrally) is violated. In order to circumvent this prob¬ 
lem we introduce the quantity 

A(t) = ln[det(i:) + tiV)] (7) 

with the pseudotime parameter t G [0,1]. For a suf¬ 
ficiently small t the approximation of Eq. ®> becomes 
valid. This property can be used together with a few 
mathematical manipulations (for details see Appendix[A|) 
to obtain the formula 


with f denoting the adjoint. 

The diagonal and the trace of an implicit matrix can 
be obtained by exploiting common probing routines |10b 
ng. A stochastic estimate of the diagonal of the linear 
operator A is given by 


1 

diag(A) = (^*~ (3) 

where * denotes a componentwise product, M = |{5}| 
the sample size, and the arithmetic mean over ^ 

with M —>■ oo. The probing vectors ^ G C" are random 
variables, whose components x (x') fulfill the condition 
= ^xx'- Analogously to the diagonal of an 
operator its trace can be probed by, e.g., 

tr(A) = . (4) 

Recently, there have been investigations to improve 
these straightforward probing methods by exploiting 
Bayesian inference m- This has been achieved by re¬ 
formulating the process of stochastic probing of an oper¬ 
ator’s diagonal (trace) as a signal inference problem. As a 
result, it requires fewer probes than the purely stochastic 
methods and thus can decrease the computational costs. 
With the phrase operator probing, be it trace or diagonal 
probing, we subsequently refer to the entirety of probing 
methods in general. 

The linear operator A can be split into a diagonal ma¬ 
trix D G C"’^” and a matrix N G which contains 

the off-diagonal part of A, i.e., 

A = DAN. (5) 


We are now interested in the value of its determinant or of 
its log-determinant, A = ln[det(A)]. In case A is mainly 
dominated by its diagonal (i.e. ND~^ ^ 1 spectrally), 
a Taylor expansion of the log-determinant might be a 
reasonable approximation. 


A = ln[det(D -b N)] 

= ln[det(D)] -b tr [ND-^] + O (tr {ND-^Y ) , 


( 6 ) 


which is sometimes feasible dealing with implicit oper¬ 
ators, e.g., see Refs. for recent applications in 

cosmic microwave background physics. This approxima¬ 
tion, however, breaks down when the relation ND~^ <C 1 


A = / dt tr 


N{D + tNy 


■A(0) 


J\t (et7v(i? + t7V)-'e)^^j+A(0) 


( 8 ) 


that represents a stochastic estimate of the log- 
determinant of A using operator probing. In particular, 
the following steps are required to evaluate Eq. ([^; 

1. Diagonal (operator-) probing to split A into 


A = diag(A) -b A — diag(A), 

= D =N 


2. an approach to invert D + tN in Eq. e.g., the 
conjugate gradient method m, 

3. trace (operator-) probing to evaluate the integrand, 

4. a numerical integration method, e.g., applying 
Simpson’s rule. 


It might immediately strike the eye of the reader that 
one recaptures the simple first-order Taylor-expanded 
version of the log-determinant, Eq. when dropping 
the pseudotime dependency of the integrand in Eq. ([^ 
by requesting t = 0. This means that in case of deal¬ 
ing with diagonal dominant operators the value of the 
correct log-determinant might be received by a coarse 
numerical integration since the integrand close to t = 0 
already yields the main correction, which might decrease 


the computational costs, see Sec. |IIB 

Equation Q further represents the main result of this 
paper and can be regarded as a spec ial case of calculat¬ 
ing partition functions (see Sec. Ill and Refs. [TBl ITT]!. 
Although the first line of it, the integral representation of 
the log-determinant, was also, independently of our work, 
found by mathematicians 10 years ago [iH], it is (to our 
knowledge) not known in the community of physics or sig¬ 
nal inference. The connection to stochastic estimators, 
however, is a novel way to evaluate the log-determinant 
of implicitly defined matrices that enables previously im¬ 


possible calculations, see Sec. Ill 


B. Numerical example 

We address here a simple and also exactly solvable nu¬ 
merical example referring to (Bayesian) signal inference 
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problems or, in general, statistical problems in physics 
(see Secs. Ill A and IIIBl, where the log-determinant 
of a covariance matrix A is of interest. If we assume 
statistical isotropy and homogeneity of a physical field, 
its covariance matrix can be parametrized by a so-called 
power spectrum. This is often a reasonable assumption^, 
e.g., in astronomy and physical cosmology, when apply¬ 
ing the cosmological principle. In this case the covariance 
matrix becomes diagonal in Fourier space. 


^kk' — ^k^kk' 


(9) 


with respective Fourier modes fc, k' and power spectrum 
Cfc. It is straightforward to show that the position space 
representation of Akk'^ given by = I’l^Akk'^k'x' 

with Fourier transformation is nondiagonal if and only 
if Cfc 7 ^ const yk. In order to apply the stochastic esti¬ 
mator of the log-determinant we use two special forms of 
the power spectrum, given by 

^ Jl+W 

with a set to 2 or 4. A value of a = 2 describes 
a mostly diagonal dominant matrix, whereas a = 4 
exhibits a significant nondiagonal structure in position 
space. To be precise, in the following we use a regular, 
two-dimensional, real-valued grid (overT^) of n = 20 x 20 
pixels to represent our position space, resulting in a ma¬ 
trix A consisting of n x n = 1.6 x 10® real numbers. See 
Fig.[T] for an illustration thereof. 

For both matrices, which we refer to as A 2 and A 4 , 
we apply Eq. ([^ given an explicit and implicit numerical 
implementation. For the explicit variant there also exist 
well-understood, precise numerical methods^ to calculate 
the determinant. Therefore, the numerical results of such 
a method can be regarded as our gold standard and hence 
serve as a reference for the probing results. Henceforth 
we will refer to it using the subscript “correct”. Both 
variants, the explicit and implicit implementation, are 
realized using the tools of NIFTy [T^ . 

After the separation of A 2 and A 4 into diagonal and 
off-diagonal parts by applying diagonal probing we calcu¬ 
late the integrands of Eq. <§ for the TO-part-discretized 
interval of t S [ 0 , 1 ] by using the conjugate gradient 
method as well as trace probing and perform the numeri¬ 
cal integration afterwards by using Simpson’s rule. The 
operator probing as well as the conjugate gradient 
method have also been realized using NIFTy. Further¬ 
more we introduce the quantities 

A(x) = / dt iv N[D + tN)~^ -|-A(0), a; G [0,1] 
do 

( 11 ) 


^ Referring to Bayesian evidence calculations such a matrix might 
be the prior or posterior covariance, see Sec. Illllfor details. 

^ See, for instance, the method described at http: 
//docs.scipy.org/doc/numpy/reference/generated/numpy. 
linalg.slogdet.html which is based on LU-factorization. 



- 0.2 4.7 


FIG. 1: Illustration of the matrices A 2 (top) and ^4 (bottom) 
in position space with linear color bars. 


to study the convergence to the final value and A(m) to 
investigate the dependency on the discretization of the 
integration interval, see Figs. andjg 


We used a rather low sample size of M = 8 for trace 
and diagonal probing [see Eqs. 0 and ([^] to demon¬ 
strate the applicability of the method to large data sets. 
The discretization of the pseudotime interval into m parts 
was chosen to be m = 10 ® for A 4 and only m = 10 for A 2 , 
see in particular Fig. [^ which illustrates the dependence 
of the probing result on m. 
















4 


TABLE I: Results of the numerical determinant calculations 
with and without probing. The absolute errors of the probing 
method are defined by ei = |Aexpiicit(l) - Aimpiicit(l)| and 
£2 = lAcorrect “ Ainipiicit(1)|• Differences between ei and £2 
arise from the discretized, numerical integration. 



^2 

A 4 

A(0) 

-1308.05 

-1771.57 

^correct 

-1566.99 

-3107.28 

^explicit(l) 

-1566.81 

-3107.29 

^implicit(l) 

-1565.33 

-3108.41 

m 

10 

1000 

M 

8 

8 

ei 

1.48 

1.12 

£2 

1.66 

1.13 


C. Discussion 


The exact numerical values of the determinant calcu¬ 
lation using explicit and implicit representations of A 4 
and A 2 can be found in Table |lj The results of the prob¬ 
ing method (implicit) compared to the correct and the 
explicit method, where Eq. Q can be evaluated without 
using a conjugate gradient or probing techniques, are 
accurate for both matrices. It is remarkable that despite 
using a relatively small sample size of M = 8 the abso¬ 
lute errors remain relatively small. The reason for this 
is that the pseudotime integration over all probed inte¬ 
grands averages the probing error. This is of particular 
importance when applying the log-determinant probing 
to large data sets, where a large sampling size should be 
avoided to save computational time. These errors can be 
decreased further, of course, by an increase of the sam¬ 
pling size and a refinement of the numerical integration. 

The results of the trace (integrand) probing and the 
determinant’s convergence behavior as well as their re¬ 
spective errors with respect to the explicit representation 
can be found in Figs.j^an di Note that the scaling of the 
ordinate is logarithmic. For both matrices, but especially 
for A 4 , the largest contribution to the integral of Eq. ^ 
comes from late t-values. Therefore, if dealing with big 
data sets, one could divide the integration interval not 
into m equal parts but by starting with a rather coarse 
discretization for small t-values and subsequently refin¬ 
ing it for larger values, e.g., by substituting dt by d\n{t') 
and thereby saving computational costs. This, however, 
might depend on the particular shape of the matrix and 
has to be studied case by case. 

The dependency of the numerical value of the deter¬ 
minant of A 4 on the discretization (in m equal parts) 
of the integration interval can be found in Fig. and 
shows that even a small value of m adds significant cor¬ 
rections to the result. The result for m oc 0(10) is, for 
instance, better than just using the determinant of the 
diagonal, A(0). This might be used in practice to in¬ 
vestigate cheaply whether the nondiagonal structure of a 
matrix influences the determinant significantly. 



0.0 0.2 0.4 0.6 0.8 1.0 

pseudo time t 


FIG. 2: The integrand of Eq.|^(lower panel) and A(a;) (upper 
panel) for explicit and implicit representations of A4. 


A huge advantage of the probing method discussed 
here is the possibility to parallelize the numerical calcula¬ 
tion almost completely. To be precise, the diagonal prob¬ 
ing beforehand, the pseudotime integral, as well as every 
single trace probing can be parallelized fully. The only 
operation that cannot be parallelized is the conjugate 
gradient method as it is a potential minimizer, using at 
least the previous step to calculate the next one. 

The determination of a suitable choice of the involved 
parameters m and M as well as the precision parameters 
for the used conjugate gradient approach and numer¬ 
ical integration method depend highly on the matrix to 
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pseudo time t 


FIG. 3: The integrand of Eq.|^(lower panel) and A(a;) (upper 
panel) for explicit and implicit representations of A 2 with only 
m = 10 steps in pseudotime. 


be studied. The computational costs and precision of 
the introduced determinant calculation thus depend on 
the combination of the chosen methods for diagonal and 
trace probing, numerical integration, the method to nu¬ 
merically invert the matrix D + tN, and the matrix A 
itself. Since it is therefore not possible to make general 
statements we consciously avoid here such a discussion 
of computational costs and precision with respect to m 
and M. A more pragmatic way to estimate these param¬ 
eters would be to downscale the problem of interest until 
the matrix of interest fits into the memory of the com- 



FIG. 4: Dependency of the determinant’s result on the dis¬ 
cretization of the integration interval into m parts, using A4. 

puter and to subsequently perform mock tests to obtain 
a suitable choice for the parameters discussed above. Af¬ 
terwards these values can be extrapolated to the size of 
the real problem. 

III. APPLICATIONS IN SCIENCE 

Within this section we present a selection of possible 
applications in science. Although there are a vast num¬ 
ber of research fields and topics which might benefit from 
the stochastic estimation of a log-determinant we focus 
henceforth on a selection of usages in Bayesian signal 
inference, in particular in physics and only present sim¬ 
ple examples. Exact, more complicated examples can be 
found in the cited works within this section. 


A. Evidence calculations & model selection 

The Bayesian evidence V{d) is a measure for the qual¬ 
ity of the model and hence for all assumed model param¬ 
eters for the data d [20]. To keep it short and simple 
we assume a model that describes a linear measurement 
of a Gaussian signal s with additive, signal-independent, 
Gaussian noise n, i.e., 

d = Rs + n, (12) 

where R represents a linear operator. A Gaussian distri¬ 
bution of a variable x is defined by 

V{x) = g{x,X) = (13) 

with related covariance matrix X and mean 

{^)v{x) ^ j'Dxx V{x). (14) 
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f V[\ denotes a phase space integral and | • | the deter¬ 
minant. Under these circumstances the evidence can be 
calculated as 


V{d) = j Vs j Vn V{d, 


t) 


^ J Vs JvnMd-Rs-n)PHs)PU) 


' \27TC,\d\ 

|27rC's||27rC'„ 


exp{-i j1'Cs|dj)|, 

(15) 


with 


J = R^C-^d, 




(16) 


and the signal and noise covariances Cs and C„, respec¬ 
tively. Therefore, to calculate the Bayesian model evi¬ 
dence, one often^ has to calculate determinants of covari¬ 
ance matrices. This might be done by probing [Eq. m 
if dealing with implicit matrices [last line of Eq. (15)] 
instead of performing the multidimensional integral sec¬ 
ond last line in Eq. s numerically. The latter has 
been done, for instance, in the field of inflationary cos¬ 
mology by the method of nested sampling |231l24j . 

This is especially of importance in the field of model 
selection or comparison | 20 | , where from an observation - 
the data - one wants to infer which theory reproduces the 
observation best. Switching from one model to another 
means, for instance^, to exchange R in Eq. (15), which 


directly affects the determinant containing Cs\d- Thus, 
the calculation of the determinant is mandatory here. 


procedure to circumvent this problem is to marginalize 
over these parameters, 

'P{jp^\d) (X J Vs J Vn V{d,s,n\p’‘)V{p^). (17) 

To continue with the simple example of Sec. ImA] we 
assume again Gaussian distributions for s and n and a 
linear measurement but with explicit dependency on p*, 
i.e., d = (i?s)[p*] -I- n. If we further follow the example 
of calibration, the parameter p* might be a calibration 
coefficient, thus affecting only R. This yields (i?s)[p®] = 
i?[p®]s and therefore 

VW) xi^jvsg{d- R[p^]s, c„) 5(s, C«)| V{p^). 

( 18 ) 

This integration can be performed analytically, produc¬ 
ing an in general non-Gaussian probability distribution 
with p®-dependent normalization (and exponent) similar 
to Eq. (15), 


V[p''\d) oc Y^|27rC'^|d[p*]| 'P(p*) 


X exp<j 


(19) 


with Cs^dW] and j[p®] now containing i?[p®] instead of R. 
In case the covariance matrices or i?[p*] are only given by 
a compnter routine (implicit representation of a matrix) 
one could use Eq. ^ to probe the determinant. 

A variety of scientific fields are affected by this prob¬ 
lem. For example, the extraction of the level of non- 
Gaussianity of the cosmic microwave background UlM 
in cosmology, the problem of self-calibration [26H28] in 
general, or leasing in astronomy |29j . 


B. Posterior distribution including 
marginalizations 


C. Realistic astronomical example 


In the field of signal inference one is typically inter¬ 
ested in reconstructing a set of i parameters p* with un¬ 
certainty from some observation, the data d. This in¬ 
formation is delivered by the posterior, given by® |25] 
7^(p*|d) cx 7^(d|p*)7^(p*). Often, however, this inference 
problem is degenerate, caused by a so-called nuisance pa¬ 
rameter. For example, consider the calibration of an in¬ 
strument is of interest and not the signal. In this case the 
signal s represents the nuisance parameter. The common 


® By the word “often” we refer to cases, in which at least one 
marginalization [see Eq. l |17| l] can be performed analytically (ap¬ 
proximated with high precision) to obtain a model-dependent 
determinant. 

^ We focus here on R for simplicity only. One could also, addi¬ 
tionally, exchange the prior covariances Cn and Cs, the assumed 
prior statistics, the parametrization of the data, and so on. 

® Note that in this case the evidence is just a scalar which normal¬ 
izes the posterior, therefore we merely state proportionalities. 


In order to study a more realistic example we consider 
a measurement device with spatially constant but un¬ 
known calibration amplitude, parametrized by I-I -7 G R, 
scanning a specific patch of the sky. The measured and 
assumed to be Gaussian sky signal s is affected by the 
instrument via a convolution C with a Gaussian kernel of 
standard deviation a = 0.05. Additionally, the observa¬ 
tion might be disturbed by fore- and backgrounds. For 
this reason we include an observational mask Mq, which 
cuts out 20% of the sky. The noise n is still assumed 
to be Gaussian and uncorrelated with the signal. Hence, 
the measurement equation is given by 

d = RIj] s + n = (1 -I- 'y)MoCs + n. ( 20 ) 

To calibrate the measurement device the calibration pos¬ 
terior V{'y\d) has to be determined. The resulting cal¬ 
ibration mean {'j)-p(^^\d) t)e regarded as an external 
calibration if the a priori knowledge on the signal is suf¬ 
ficiently strong. Otherwise one could infer the signal and 
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calibration amplitude 7 simultaneously from data using 
iterative approaches [27] . Using Eq. (191 as well as a flat 
prior on 7 we obtain 

lnlP( 7 |d) = - ^ In C'q] [t] + W Cqd [7] J H ^21) 
+ const., 

which exhibits in particular the 7 -dependent determinant 


^7] [7] = 1(1 + 7 ) + 7) + c; 


For the nume rical evaluation of Eq. 
settings of Sec. 


IIB 


( 22 ) 
we use the 


with Ca{k,k') = (1-I- fc) ^ 5kk 


calibration amplitude parameter of 7 = 2 , and a noise 
covariance of [Cn)^ to generate a data re¬ 

alization. The pseudotime interval has been discretized 
into 10^ parts. The numerically determined calibration 
posterior for a given data realization can be found in 
Fig.§ which demonstrates again the efficiency of the 
stochastic method using only eight probes for a single 
trace probing operation. The figure also illustrates the 
impact of the determinant on the log-posterior, which 
would not peak in the shown interval without it. 



FIG. 5: Logarithmic posterior of the calibration amplitude 
parameter 7 using implicit and explicit representations of the 
involved operators, see Eq. (21 1 and Eq. (22 I for details. The 
abbreviation A denotes the logarithm of the term given by 
Eq. (|22|. 


IV. SUMMARY 

Motivated by the problem of finding a way to efficiently 
determine the determinant of an implicitly defined ma¬ 
trix or operator, we derived a formula, Eq. (|^, repre¬ 
senting a stochastic estimate of its log-determinant. This 
has been achieved by reformulating the log-determinant 
as an integral representation and transforming the in¬ 
volved terms into stochastic expressions, which includes 


a numerical integration and a trace probing. Numerical 
examples have shown that the discretization of the inte¬ 
gration interval may be very coarse in case the probed 
operator is sufficiently diagonal. In case it exhibits a sig¬ 
nificant nondiagonal structure one has to fine-grain the 
discretization of this interval. The number of probes nec¬ 
essary for the trace probing, however, remains very low 
in the studied examples. These facts combined with the 
almost complete parallelizability of this approach might 
keep the computational costs within reasonable limits in 
many situations. 

This method clearly has more general applications but 
might in particular be useful for Bayesian signal inference 
and model comparison when dealing with large data sets 
as often given, for instance, in astronomy and cosmology. 
To be precise, it might be beneficial in all fields where the 
numerical calculation of a determinant of an operator is 
mandatory. 
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Appendix A: Integral representation of the 
log-determinant of a matrix 

Here Eq. (|^ is derived. Following Sec. [IT] the log- 
determinant A of an operator A can be parametrized 
by A = ln[det(I? -I- N)] with D being the diagonal and 
N the off-diagonal part of A. Since A can be Taylor- 
expanded for small N (spectrally compared to D) only, 
we employ a method from the field of renormalization 
theory |28l [30] . Accordingly, we introduce an expansion 
parameter 6 t 1 to suppress the influence of N. In par¬ 
ticular, we replace A by ln[det(Z/ + i5fA)] for a moment. 
For sufficiently small values of 5t, in the following inter¬ 
preted as tiny pseudotime steps, we can approximate A 
by Eq. (§. Theoretically, a single pseudotime step could 
be infinitesimal small, enabling the formal definition of 
the derivative 

dA{t) _ ln[det(i/ -f (f -f St)N)] — ln[det(D -|- tN)] 
dt it-fO 5t 

= lim — In [det (l -|- 5tN[D 
5t—^0 ot 

= lim ^tr [in (l -|- StN[D + t7V]“^)] 

6 t —>^0 ot 

= tr [V[D + tV]-7, 

(Al) 

with the definition 

A{t)=\n[det{D + tN)]. (A2) 






















Integrating the pseudotime derivative of A(t) yields the 
integral representation of the log-determinant, 


A = 


dt tr 


N{D + tNy 


A(0). 


(A3) 


the inverse matrix of the integrand of Eq. (A31. 


Finally, we replace the trace by stochastic trace prob¬ 
ing and perform the pseudotime integral by an numeric 
integration method. This yields 


This integral representation has also been found by 
Ref. |18| . where its validity has been proven for weak 
diagonal dominant and Hermitian positive definite ma¬ 
trices. In particular one has to ensure the existence of 


A = J^dt + (A4) 
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